#ifndef PPSPLINE_H_
#define PPSPLINE_H_
 
#include<vector>
#include"polynomial.h"
#include"Function.h"
#include"linearEqnSolver.h"
#include"dividedDiff.h"

using namespace std;

enum boundary_condition{
    complete,
    second,
    natural,
};

/**
 * pp样条需要存储一组节点, 以及这组节点所形成区间上的多项式
 */
template<int order>
class ppSpline
{
private:
    vector<double> knots;
    vector<polynomial<order>> polys;
public:
    ppSpline(){};
    ~ppSpline(){};
    ppSpline(const vector<double>& _knots, vector<polynomial<order>>& _polys):
        knots(_knots), polys(_polys){};

    vector<double> get_knots() const { return knots; }

    vector<polynomial<order>> get_polys() const { return polys; }

    int get_num_of_knots() const{ return knots.size(); }    // 得到节点个数

    void print_ppForm_data(int sample_num, ostream &out);   // 打印密集散点以作图
                                                            // 传入小段区间上采样点个数

    void print_ppForm();    // 打印样条多项式

    void print_plot_data(int sample_num, string &name); // 给出文件名以存储密集散点

};

// 成员函数实现

template<int order>
void ppSpline<order>::print_ppForm_data(int sample_num, ostream &out){
    int N = get_num_of_knots();
    double delta; 
    for(int i = 0; i < N-1; i++){
        double x = knots[i];
        delta = (knots[i+1] - knots[i]) / sample_num;
        for(int j = 0; j < sample_num; j++){
            out << x << " " << polys[i].get_value_at_x(x) << endl;
            x += delta;
        }
    }
    out << knots[N-1] << " " << polys[N-2].get_value_at_x(knots[N-1]) << endl;
    return;
}

template<int order>
void ppSpline<order>::print_ppForm(){
    int N = get_num_of_knots();
    for(int i = 0; i < N-1; i++){
        cout << "区间[" << knots[i] << "," << knots[i+1] << "]:";
        polys[i].print_poly();
    }
} 

template<int order>
void ppSpline<order>::print_plot_data(int sample_num, string &name){
    ofstream file(name);
    if(file.is_open()){
        print_ppForm_data(sample_num, file);
        file.close();
        cout << "File " << name;
        cout << " created and written successfully.\n";
    }
    else{  
        cout << "Error creating file: " << name << ".\n";  
    }
}

// 给出插值点,函数,边界条件,三次样条拟合函数y=f(x)
ppSpline<3> solver( const vector<double> &knots, const Function &fun,
                    const boundary_condition boundary );

// 给出一维散点,三次样条拟合,边界条件为自然样条
ppSpline<3> solver( const vector<double> &knots,
                    const vector<double> &vals,
                    const boundary_condition boundary );

// 给出一维散点,三次样条拟合,边界条件为complete或second
ppSpline<3> solver( const vector<double> &knots,
                    const vector<double> &vals,
                    const boundary_condition boundary,
                    const double &left_bound, const double &right_bound );

// 给出高维散点,三次样条拟合,边界条件为自然样条
vector<ppSpline<3>> solver( const vector<vector<double>> &points, 
                            const boundary_condition boundary );

// 给出高维散点,三次样条拟合,边界条件为complete或second
vector<ppSpline<3>> solver( const vector<vector<double>> &points, 
                            const boundary_condition boundary,
                            const vector<double> &left_bound, 
                            const vector<double> &right_bound );

// 给出插值点,函数,一次样条(分段线性插值)拟合y=f(x)
ppSpline<1> solver(const vector<double> &knots, const Function &fun);

// 给出一维散点,一次样条(分段线性插值)拟合
ppSpline<1> solver(const vector<double> &knots, const vector<double> &vals);

// 给出高维散点,一次样条(分段线性插值)拟合
vector<ppSpline<1>> solver(const vector<vector<double>> &points);

// 打印对高维散点拟合后的密集点用以作图, sample_num为每段小区间上均匀采样的样本点数
template<int order>
void print_dis_pt_data( const vector<ppSpline<order>> &n_dim_ppSpline, 
                        int sample_num, ostream &out)
{
    int n = n_dim_ppSpline.size();
    int N = n_dim_ppSpline[0].get_num_of_knots();
    vector<double> knots = n_dim_ppSpline[0].get_knots();
    for(int i = 0; i < N-1; i++){
        double delta = (knots[i+1] - knots[i]) / sample_num;
        for(double t = knots[i]; t < knots[i+1]; t += delta){
            for(int j = 0; j < n; j++){
                out << n_dim_ppSpline[j].get_polys()[i].
                                get_value_at_x(t) << " ";
            }
            out << endl;
        }
    }
    for(int j = 0; j < n; j++){
        out << n_dim_ppSpline[j].get_polys()[N-2].
                get_value_at_x(knots[N-1]) << " ";
    }
    out << endl;
    return;
}

#endif